#ifndef _BIOM_
#define _BIOM_ 1

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <algorithm>

namespace biom
{
    template <class T>
    class matrix
    {
    public:
        std::vector<T> block;
        size_t row, column;
        std::string name;

        // Constructor
        matrix();
        matrix(const size_t r, const size_t c, const std::string &n = "");
        matrix(const size_t r, const size_t c, const std::initializer_list<T> ilist, const std::string &n = "");
        ~matrix()
        {
            row = 0, column = 0, name = "";
            block.clear();
        }

        matrix<T> &clear();

        T &operator()(size_t i, size_t j);
        matrix<T> &operator=(const matrix<T> &OBJ);
        matrix<T> &operator=(const std::initializer_list<T> ilist);
        matrix<T> operator+(const matrix<T> &B);
        matrix<T> operator-(const matrix<T> &B);
        matrix<T> operator*(const matrix<T> &B);
        bool operator==(const matrix<T> &B);
        bool operator!=(const matrix<T> &B);

        friend std::ostream &operator<<(std::ostream &output, const matrix<T> &OBJ)
        {
            output << OBJ.name << " <" << OBJ.row << ',' << OBJ.column << '>' << std::endl;
            for (size_t i = 0; i < OBJ.row; i++)
            {
                for (size_t j = 0; j < OBJ.column; j++)
                {
                    if (j != 0)
                        output << setiosflags(std::ios::right) << std::setw(8) << OBJ.block[(i)*OBJ.column + j];
                    else
                        output << '[' << setiosflags(std::ios::right) << std::setw(8) << OBJ.block[(i)*OBJ.column + j];
                }
                output << "\t]\n";
            }
            return output;
        }
    };

    template <class T>
    matrix<T>::matrix() : row(0), column(0), name("") { block.clear(); }

    template <class T>
    matrix<T>::matrix(const size_t r, const size_t c,
                      const std::string &n)
        : row(r), column(c), name(n) { block.clear(); }

    template <class T>
    matrix<T>::matrix(const size_t r, const size_t c,
                      const std::initializer_list<T> ilist,
                      const std::string &n)
        : row(r), column(c), name(n)
    {
        block.clear();
        block.assign(ilist.begin(), ilist.end());
    }

    template <class T>
    T &matrix<T>::operator()(size_t i, size_t j)
    {
        return block[j + column * i];
    }

    template <class T>
    matrix<T> &matrix<T>::clear()
    {
        block.clear();
        row = 0, column = 0, name = "";
        return *this;
    }

    template <class T>
    matrix<T> &matrix<T>::operator=(const matrix<T> &OBJ)
    {
        block.clear();
        block.assign(OBJ.block.begin(), OBJ.block.end());
        this->row = OBJ.row, this->column = OBJ.column, this->name = OBJ.name;
        return *this;
    }

    template <class T>
    matrix<T> &matrix<T>::operator=(const std::initializer_list<T> ilist)
    {
        block.clear();
        block.assign(ilist.begin(), ilist.end());
        return *this;
    }

    template <class T>
    matrix<T> matrix<T>::operator+(const matrix<T> &B)
    {
        matrix<T> RES(this->row, B.column, "(" + this->name + "+" + B.name + ")");
        for (size_t i = 0; i < row * column; i++)
        {
            RES.block.push_back(block[i] + B.block[i]);
        }
        return RES;
    }

    template <class T>
    matrix<T> matrix<T>::operator-(const matrix<T> &B)
    {
        matrix<T> RES(this->row, B.column, "(" + this->name + "-" + B.name + ")");
        for (size_t i = 0; i < row * column; i++)
        {
            RES.block.push_back(block[i] - B.block[i]);
        }
        return RES;
    }

    template <class T>
    matrix<T> matrix<T>::operator*(const matrix<T> &B)
    {
        matrix<T> RES(this->row, B.column, "(" + this->name + "*" + B.name + ")");
        for (size_t i = 0; i < this->row; i++)
        {
            for (size_t j = 0; j < B.column; j++)
            {
                T temp = 0;
                for (size_t c = 0; c < this->column; c++)
                {
                    temp += (this->block[(i) * (this->column) + c] *
                             B.block[(c)*B.column + j]);
                }
                RES.block.push_back(temp);
            }
        }
        return RES;
    }

    template <class T>
    bool matrix<T>::operator==(const matrix<T> &B)
    {
        if (row == B.row && column == B.column)
        {
            for (size_t i = 0; i < row * column; i++)
            {
                if (this->block[i] - B.block[i] > 0 || this->block[i] - B.block[i] < 0)
                    return false;
                else
                    continue;
            }
            return true;
        }
        else
        {
            return false;
        }
    }

    template <class T>
    bool matrix<T>::operator!=(const matrix<T> &B)
    {
        bool truth = this->operator==(B);
        return !truth;
    }

    template <class T>
    matrix<T> tran(const matrix<T> &Mt)
    {
        matrix<T> RES(Mt.column, Mt.row, "(~" + Mt.name + ")");
        for (size_t i = 0; i < RES.row; i++)
        {
            for (size_t j = 0; j < RES.column; j++)
            {
                RES.block.push_back(Mt.block[(j)*Mt.column + i]);
            }
        }
        return RES;
    }

    template <class T>
    matrix<T> operator~(const matrix<T> &Mt)
    {
        return tran(Mt);
    }

}

#endif